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PREFACE 


This  publication  (interim  Report  69-2)  forms  the  second  part  of 
two  reports  dealing  with  gas  exchange  from  soil.  Included  here  are  pages 
125-138  forming  the  appendix.  We  have  chosen  to  make  two  reports  because 
of  the  anticipated  diversity  of  interest  in  them. 


E.  R.  Lemon 


TABLE  OF  CONTENTS 


i  I 


Page 

3.  Relationship  between  cross-spectral  density  functions 

and  2-D  auto- covariance  functions  .  183 

VI.  Stannary . 185 

VII.  Literature  Cited  .  I87 


1: 


LIST  OF  ILLUSTRATIONS 


Figure 

A. 2.1  Illustration  of  the  sampling  of  p  at  the  various  times, 

tj . 

A. 2. 2  Plot  of  Fourier  line  spectrum  . 

A. 2. 3  Periodogram  of  •—  2 1 Cn I  2  against  f . 

Af  1  1 

A.2.h  Plot  of  spectral  density  function  against  frequency 
A. 2. 5  The  curve  for  calculation  of  a  spectrum  by  hand  .  .  . 

A. 2. 6  A  spectrum  . 

A. 2. 7  Three  methods  for  plotting  spectra  of  air  pressure  .  . 


Page 

129 

135 

137 

138 
l4l 

145 

146 


A. 2. 8  The  au^o-covariance  function  for  the  curve  in  Fig.  A. 2. 5  .  .  .  148 

A. 2. 9  Illustration  cf  hov  digital  sampling  may  cause  a  high  frequency 

component  to  alias  and  appear  as  a  low  frequency  component  .  .  1.55 

A. 3.1  Coherence  and  phase  spectrum  for  air  pressure  between  a  point 
on  the  ground  surface  in  a  field  and  another  point  downwind 


from  the  first . l6l 

A.  4.1  Illustration  of  coordinate  system  and  notation  for  sampling 

p . 165 

A.  4.2  Plot  of  Fourier  line  spectrum . l69 


A.  4.  3  Plot  of  3-D  "periodogram" . . . 171 

A. 4.  4  Contour  surface  of  spectral  density  function, 

plotted  against  frequency  and  wave  number  .  172 


A. 5.1  Sampling  scheme  which  can  be  used  where  p  is  statistically 
stationary  and  spacially  homogeneous  . 


179 


Appendix 


BASIC  CONCEPTS  OF  SPECTRAL  ANALYSIS  BY  DIGITAL  MEANS 

I.  INTRODUCTION 

A  statistical  technique  that  has  helped  meteorologists  study  wind 
and  other  fluctuating  variables  but  which  as  yet  has  been  rarely  used  by 
agronomists  is  spectral  analysis.  Spectral  analysis  is  used  to  evaluate, 
the  contributions  of  different  frequencies  of  the  fluctuation  to  the 
total  variance  of  an  entity,  such  as  wind  velocity,  which  changes  with 
time.  A  graph  of  the  variance  per  frequency  band  forms  a  spectrum  not 
unlike  the  more  cocmionly  known  spectrum  of  light  intensity  versus  wave 
length.  Knowing  the  frequencies  of  the  dominant  modes  of  oscillation  (or 
knowing  that  there  are  none)  can  help  one  to  understand  and  visualize  the 
physics  of  the  transport  processes  mentioned  in  the  first  paragraph. 

Allen  (1968),  for  instance,  has  used  spectral  analysis  to  study  the  eddy 
structure  of  wind  in  a  Japanese  larch  canopy  and  found  significant  contri¬ 
butions  to  the  variance  of  wind  velocity  by  eddies  the  size  of  the  tree 
spacing  at  mid-canopy  heights.  There  was  surprisingly  little  contribu¬ 
tion  from  these  eddies  at  the  bottom  of  the  canopy. 

Spectral  analysis  car.  be  used  to  study  two  or  more  fluctuating  vari¬ 
ables  to  determine  the  closeness  of  their  relationship  for  different 
frequencies  of  change,  Desjardins  (1967),  for  instance,  has  compared  the 
responses  of  various  atmometers  to  open  pan  evaporation.  He  found  that  the 
variances  were  similar  for  periods  of  time  longer  than  five  days,  but 
different  for  shorter  periods,  indicating  that  different  physical  factors 
must  be  in  operation  for  the  different  instruments .  Rodrigue z-Iturbe  and 
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Norain  (1969 )  used  spectral  analysis  techniques  to  show  the  existence  of 
strong  correlations  among  annual  oscillations  in  precipitation  and  run-off 
in  the  Pacific  Coast  region  of  the  United  States.  The  annual  cycle  of 
precipitation  could  be  considered  the  same  for  observation  stations  up  to 
1000  km  apart. 

Spectral  analysis  was  first  developed  for  practical  use  by  electri¬ 
cal  engineers.  Because  the  variance  of  a  voltage  across  or  a  current 
through  a  unit  resistance  is  proportional  to  bhe  average  power  dissipation, 
the  term  power  spectrum  has  often  been  used  to  refer  to  the  spectrum 
of  variance.  Communication  engineers  have  made  much  use  of  spectral 
analysis,  and  the  theory  and  procedures  of  making  spectral  analysis  are 
explained  in  much  detail  by  Blackman  and  Tukey  (19!;  8).  A  much  less  sophis¬ 
ticated,  but  easily  understood  explanation  of  the  underlying  principles 
of  spectrum  analysis  from  the  point  of  view  of  meteorology  may  be  obtained 
from  Panofsky  and  Brier  (1965).  Luraley  and  Panofsky  (1961+)  present  a 
much  more  sophisticated  treatment  from  the  same  viewpoint.  Jenkens  and 
Watts  (1968)  present  an  excellent  explanation  of  the  principles  from  a 
statistician's  point  of  view. 

In  1965  Cooley  and  Tukey  revealed  a  new  algorithm  for  fast  computa¬ 
tion  of  Fourier  transforms  which  has  made  spectral  analysis  much  easier. 

An  easily  understood  description  of  the  new  method  is  given  by  Brigham 
and  Morrow  (1967).  More  information  may  be  obtained  in  the  June  1967 
issue  of  IEEE  Transactions  on  Audio  and  Electroacoustics,  vol.  15,  no.  2, 
which  was  devoted  entirely  to  the  fast  Fourier  transform  and  how  it  re¬ 
lates  to  spectral  analysis  and  other  subjects. 
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Hie  author  has  found  all  of  the  above  references  useful  in  under¬ 
standing  spectral  analysis,  but  found  none  which  developed  the  principles 
with  reference  to  the  new  computational  methods  for  persons  who  previously 
had  little  knowledge  about  spectrum  analysis.  This  paper  attempts  to 
develop  the  concepts  of  spectral  analysis  for  a  reader  who  previously 
knows  little  about  spectral  analysis  and  whose  background  includes  no  more 
mathematics  than  basic  calculus  and  statistics. 

The  orientation  here  will  be  toward  the  use  of  a  digital  computer  to 
compute  spectra  from  a  series  of  samples  of  discrete  data;  therefore,  the 
discrete  forms  of  equations  are  used  wherever  possible.  Spectra  can  be 
measured  directly  by  electronic  filtering  of  continuous  electrical  analog 
signals,  but  suitable  analog  sisals  cannot  be  obtained  for  many  varia¬ 


bles  . 
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II.  ONE-DIMENSIONAL  SPECTRAL  ANALYSIS 

This  section  will  discuss  the  fluctuation  of  an  entity  that  changes 
with  time.  All  of  the  concepts  would  apply  equally  well  if  the  entity 
would  change  with  some  other  variable,  such  as  distance,  but  spectral 
analysis  has  been  used  more  frequently  with  entities  that  change  with 
time  than  with  other  variables,  so  this  orientation  will  be  used  in  this 
development.  The  term,  one-dimensional,  refers  to  the  fact  that  the 
entity  fluctuates  only  with  one  variable,  time.  Later,  an  entity  that 
fluctuates  with  two  variables,  time  and  distance,  will  be  examined,  and 
the  analysis  will  be  referred  to  as  two-dimensional  spectral  analysis. 


1.  Fourier  series  representation 

Suppose  some  quantity,  p  (perhaps  air  pressure),  a  function  of  time, 
t,  has  been  observed  N  times,  where  N  is  chosen  to  be  odd  for  convenience, 
and  that  the  observations  have  been  spaced  equally  At  units  apart.  Let 
the  individual  time  of  sampling  be  labelled  t^  =  j  At  for  j  -  -  , 


-  (N-3)  o  i 

- n . . . .  u  9  J-  j 


,  so  that  the  time  origin  is 


in  the  middle  of  the  total  sampling  period,  T.  Note  that  T  =  (N-l)  At. 
Figure  A. 1.1  illustrates  the  sampling  scheme. 

Since  we  are  basically  interested  in  evaluating  the  amplitude  and 
frequency  of  the  fluctuations  in  p,  it  is  useful  to  describe  p  by  a 
Fourier  series  of  complex  exponentials  as  in  Equation  A. 2.1,  where  i  =/  -1 


p(t)  =  E 

n  =  -°° 


C  expfi  n  2  ti  t) 

n  IJI 


A. 2.1 


Equation  A. 2.1  represents  the  continuous  pressure  between  time  -T/2  and 
+T/2  as  a  function  of  continuous  time.  The  complex  coefficients  Cp  may 
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be  evaluated  by  multiplying  both  sides  of  Equation  A. 2.1  by  enpf1  3.  -  11  tj 

T 


for  q  = . -2,  -1,  0,  1,  2 . and  integrating  from  -T/2  to 

+T/2.  Noting  that 
T/2 

/  exp(i-B-gJLj)  expfi  q  2  *  t)  dt 
-T/2  T  T 


=  T  when  q  =  -n 

=  0  when  q  4  -n  A. 2. 2 

one  obtains 

T/2 

C  =  1/T  /  p(t)  n2U)  dt  A.2.3 

n  -T/2  T 


Equations  A. 2.1  and  A.2.3  are  for  continuous  variables.  To  put 
Equation  A. 2.1  in  a  form  for  handling  a  finite  number  of  discrete  observa¬ 
tions,  one  writes 


N-l 

p(t.)  =  T  C  expC1  n  2  1*1) 

j  n  =  -(N-l)  n  T 

2 


A.2.U 


Equation  A.2.3  is  discretized  by  numerical  integration  with  the 
rectangle  rule  to  give 


N-l 

2 


■l^(N-l) 

2 


p(t  )  exp(~i  n  2  n  tj) 

J  ip 


A. 2. 5 


Each  pair  of  terms  in  Equation  A. 2.1+  corresponding  to  n  and  -n  repre¬ 
sents  a  particular  harmonic  in  the  fluctuations  of  p.  Also  n  denotes  the 
frequency  of  the  harmonic  because  it  measures  the  number  of  complete 
cycles  the  n*"*1  harmonic  executes  in  time  T.  Hie  particular  limits  on  n 
arise  because  of  an  assumption  inherent  in  making  the  sampled  points  of  p 
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substitute  for  the  real,  continuous  p.  This  assumption  is  that  the  sample 
interval,  At,  is  small  enou^i  so  that  there  are  no  undetected  higji  fre¬ 
quency  fluctuations  in  p  between  the  points.  As  discussed  by  Blackman  and 
Tukey  (1958,  P*  30),  the  higiest  frequency  or  harmonic  which  can  be 
detected  when  sampling  at  discrete  intervals  is  given  by  l/(2At),  which 

corresponds  to  n  = ±  N-l  .  If  p  has  higher  frequency  fluctuations,  the 

~2~ 

data  will  be  aliased.  This  problem  is  discussed  more  fully  later.  Hie 
term  corresponding  to  n  =  0  is  the  mean  of  p.  Hie  terms  corresponding  to 
n  =  ±1  represent  the  first  (or  fundamental)  harmonic  which  has  a  frequency 
of  1/T.  It  is  the  lowest  frequency  of  fluctuation  in  p  which  can  be 
distinguished,  and  if  lower  frequencies  are  present,  such  as  is  the  case 
when  there  is  a  trend  in  p,  corrections  must  be  made,  as  will  be  discussed 
later . 

Note  that,  in  general,  the  Cn  are  complex.  If  one  writes  Equation 
A.2.h  in  the  form 


N-l 

2 


p(t  )  =  C  +  2  1°  expf1  n  2  71  t3)  +  C  expC-1  n  2  17  tJ  ]  A.2.6 

.1  .  _  -i  n  m  n 


n=l 


and  uses  the  Euler  equations 


e^0  =  cose  +  i  sine 


e-^9  =  cos6  -  i  sine 


Equation  A. 2. k  is  changed  to  a  Fourier  series  of  sines  and  cosines. 


A. 2. 7 
A. 2. 8 


N-l 

2 


p(t.)  -  c  +  2 


fn  2  it 


O  -  Kcn  +  C-n)  cos^ - 

n=l  T 


iL) 


i(Cn  -  C_n)  sin(n  2  v  t ^ ) ] 
T 


A.2.9 
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Then,  if  one  defines 


So  =  2Cc 


A.2.10 


=  Cn  +  C-n 


A. 2. 11 


b  =  i  (C  -  C  ) 
n  n  -n 


A. 2. 12 


Equation  A. 2. 9  becomes  the  more  familiar 


p(t.)  =  *  +  T  [a  cosiillil)  +  K  sin  Cj.  )  ] 

J  2  n=l  T  T 


A.2.13 


which  will  not  be  used  here  because  the  Fourier  transforms  discussed  later 


in  the  paper  are  more  easily  handled  by  the  series  of  exponentials. 
Also,  note  that  by  rearranging  Equations  A. 2. 11  and  A. 2. 12, 


Cq  =  fn  -  i  ^n 

2  2 


C_n  =  %.  +  i 


A.2.10 


A. 2. 15 


so  that  Cn  and  C_n  are  seen  to  be  complex  conjugates,  a  result  which  is 


useful  later. 


2.  Relation  of  variance  to  spectral  density  function 

A  common  statistical  parameter  used  to  measure  variability  is  the 
P 

variance ,  s  .  It  is  defined 


2  _  1 
s  "  FT 


N-l 

i  !?(*,)  -  pi£ 

J  -  2 


A. 2. 16 


where  the  p(tj)  are  the  observations  of  p,  N  is  the  number  of  observations 


which 


and  p  is  the  mean  of  p(t.).  The  limits  of  summation  are  ±(N~1) t 

■J  2 

will  be  the  limits  of  summation  from  the  point  on,  unless  noted  otherwise. 

If  Equation  A. 2. 4  is  used  to  substitute  for  p(tj)  in  Equation  A.2.16, the 

variance  becomes  explicitly  expressed  as  a  function  of  the  harmonics  of  p. 

Remembering  that  p  =  C0  and  -  ±.  =  ,  one  obtains 

N-l  T 


=  ~  I  [  z  [<*  exp(LliJL!l)]  -  Co f 

T  J  n  T 


A.2.17 


Equation  A.2.17  is  the  numerical  form  for  the  following  integral  equation. 


=  1/T  /  [I  [C  exp(UL2JLt.)]  -0o]2  dt 
-T/2  n  n  T 


A.2.18 


If  the  summation  is  written  out  in  a  series,  and  if  the  squaring  is 
performed  on  the  series,  the  following  expression  is  obtained 


T/2 

’  1/T  I  2 '  +  sc-(K-i)c-(»-3)^(--(1,-f-  -*> 


+ . 2C-(N-1)  CH-1  + . 

2  2 

+  cfn  exp(zi-iLJiJLl)  +  ?C_riC_(ii+l)e*pf-i(2n  +  l)  2  *  t) 


+ . +  2C_nCn  + . +  c£  exp(i^LJLl_£.) 

+ . 

+ . ]  dt 

A. 2. 19 


Recalling  Equation  A. 2. 2,  the  variance  is  seen  to  be  equal  to  a  sum 


of  squared  moduli  of  coefficients. 


2 


13*t 


s2  = 


N-l 

2 

y  c  c 

n=l  "nn 


N-l 

2 


A. 2. 20 


Hie  vertical  bars  denote  the  modulus,  which  is  obtained  by  multi¬ 
plying  CQ  by  its  complex  conjugate,  C_Q,  and  taking  the  square  root  of 

the  product.  In  Equation  A. 2. 20  the  variance  is  expressed  as  a  function 

N-l 

of  the  amplitudes  of  the  harmonics  from  1  to  .  The  variance  may  also 
be  considered  as  a  sum  over  negative  and  positive  harmonics  since  the 
modulus  is  an  even  function.  Thus, 

N-l 


2  2  2 

s  =  -c  +  y 

o  ■*— 


A. 2. 21 


n  =  - 


(N-l) 


The  form  in  Equation  A. 2. 20  is  most  common,  however. 

Now,  if  one  defines 

f  =  n/T  A. 2. 22 

n 

and  Af  =  l/T  A.2.23 

and  plots  2jCn|2  against  fn  we  obtain  the  Fourier  line  spectrum  shown 
in  Figure  A. 2. 2.  The  f  represent  the  frequencies  of  the  harmonics,  and 
Af  represents  the  frequency  increment  between  successive  harmonics.  Thus, 
the  "curve"  is  a  series  of  lines  of  height  2|cJ2  spaced  Af  units  apart. 

If  one  assumes  that  heights  of  2 | Cn | ^  are  distributed  uniformly  over 
the  frequency  band  from  fn  -  Af  to  fn,  a  histogram  may  be  obtained 
where  the  height  of  the  bands  is  determined  from 


2 [  Cn |  2  =  (heigit)  (  Af) 


A. 2.21* 


or  height  =  —  2 1 C  | 2  A. 2. 25 

Af 

A  plot  of  —  2 1 C  { 2  against  f  is  shown  in  Figure  A. 2. 3. 

Af  n  n 

Figure  A. 2. 3.  is  called  a  periodogram  in  spite  of  the  fact  that  the 
axis  is  frequency  and  not  time.  The  total  area  of  the  rectangles  equals 
the  total  variance.  As  the  period,  T,  increases,  Af  decreases  ay  an 
inverse  proportion,  and  the  rectangles  in  the  histogram  of  Figure  A. 3 
become  narrower.  In  the  limit  as  T  -*■  00 ,  Af  tends  to  zero,  and  the  rectangles 
become  so  dense  that  their  discrete  upper  edges  approach  a  smooth  curve. 

The  limiting  smooth  curve  is  the  spectral  density  function,  s(f),  and 
the  plot  of  it  against  frequency  in  Figure  A.2.U  represents  a  variance 
spectrum  or  power  spectrum.  In  practice,  of  course,  one  does  not  have  an 
infinite  sampling  time,  as  is  implied  by  letting  T  -*■  »,  so  one  must  approx¬ 
imate  the  spectral  density  function,  using  finite  sampling  times  and 
finite  Af's. 

The  practical  spectral  density  function,  now  defined  as 

s(f  )=-i2|c  |2  n-1,  2, . 5=i  a.2.26 

Af  2 

represents  the  variance  in  p  per  frequency  increment.  The  total  area 

under  the  curve  represents  the  total  variance.  The  graph  can  be  normalized, 

if  desired,  by  divid-’'ng  s(f  )  by  the  total  variance  calculated  directly 

n 

from  Equation  A.2.1o.  The  normalized  graph  will  always  have  the  area 
under  the  curve  equal  to  1.0. 


3.  Relation  of  spectral  density  function  to  Fourier  transform 
Using  Equations  A. 2. 22  and  A. 2. 23,  Equations  A.2.2U  and  A.2.25  may 


be  rewritten  in  the  forms : 

p(t.)  =  Z  C  exp(i2Trf  t.) 
J  _  n  n  j 


A.2.27 


and 


Cn  =  AtAf  Z  p(tj)exp(-i2Trfntj) 

J 

Substituting  Equation  A.2.28  in  A.2.27,  one  obtains 


A.2.28 


p(tj)  -  Z  [AtAf  Z  p(tj  )exp(-i2irfntj  )  ]  exp( i2irfnt ^ )  A. 2. 29 

n  j 

If  one  now  defines 

N-l 

2 

S(fn)  =  At  p(tj)eXp(-i27TfntJ)  ,  A.2.30 

2 

where  the  S(fn)  are  commonly  called  Fourier  coefficients, 

N-l 

2 

then  p(tj)  =  Af  2  S(fn)exp(i2irfntj )  A.2.31 

2 

S(fn)  and  p(t^)  form  a  discrete  Fourier  transform  pair.  Using  the 

computer  technique  of  fast  Fourier  '•■rans formation ,  as  described  by  Brigham 

and  Morrow  (1967),  S(fn)  can  be  rapidly  computed  from  the  set  of  data 

points  (t,).  Alternatively,  if  the  s(f  )  are  kncwn,  the  p(t.)  can  be 
tJ  ^  J 

computed  from  the  S(fn). 

The  significance  of  the  S(f  )  may  be  realized  by  noting  that 


ll+o 

and  remembering  that  and  C_n  are  complex  conjugates.  Then  if  S(fn) 
is  multiplied  by  its  complex  conjugate,  S(-fn),  one  obtains 

S(f  ).S(-f  )  =  |S(f  )|2  =  °n  ■  c-n  =  1  |cj2  A. 2. 33 

“  “  (Af)2 

where  the  vertical  bars  denote  modulus. 

The  spectral  density  function  defined  by  Equation  A. 2. 2 6  is  seen  to 

be  equal  to  twice  the  squared  modulus  of  the  Fourier  transform  of  the  date 

points  times  the  frequency  increment. 

s(f  )  =  2  Af|s(f  )|2  ,  n  =  1,  2,  .  .  A. 2. 3b 

n  n  2 

Thus,  the  spectral  density  function  may  be  computed  by  a  fast  Fourier 

transform  of  the  data  for  n  =  1,  2  .  .  .  and  then  multiplying  each 

2 

of  the  S(fn)  by  its  complex  conjugate  times  2Af. 

An  example  of  the  use  of  Equations  A. 2. 30  and  A.2.3b  is  provided  by 

the  hand  calculation  of  the  spectrum  for  the  simple  curve  illustrated 

in  Figure  A. 2. 5.  Noting  that  At  =  0.1  sec,  that  N  =  IT,  and  that  f  t  = 

J 

nj/(N  -  l),  one  writes 

+8 

s(f_)  =  (0.1)  E  p(t.)  (cos^-ESJ-  -  i  sinfcHHiL  ) 

n  j=-8  j  16  16 

The  next  step  is  to  compute  all  the  values  of  cos  — 7Tl^-  and  sin— for 

16  16 

j  =  -8  .  .  .  ,  +8  and  n  =  1,  as  has  been  done  in  the  third  and  fourth  columns 

of  Table  A. 2.1.  In  the  fifth  and  sixth  columns,  the  cosine  and  sine  values 

are  multiplied  by  the  corresponding  p(t.)  and  summed  on  j.  The  calcu- 

J 

lation  of  S(f^)  and  s(f^)  is  illustrated  at  the  bottom  of  Table  A. 2.1. 

The  process  is  continued  for  n  =  2,  as  shown  in  the  table,  up  to  n  =  8  to 


Table  A. 2.1  Calculation  of  the  spectral  density  for  n  =  1,  2  for  the  curve  in  Pig.  A. 2. 5 


707  0.707 
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obtain  the  spectrum  illustrated  in  Fig.  A. 2. 6.  The  spectrum  shows  a 
peak  at  a  frequency  of  1.88  cycle/sec  as  one  would  predict  from  the 
distinct  repetition  near  this  frequency  in  the  original  data. 

The  spectrum  in  Figure  A. 2. 6  is  plotted  with  a  linear  scale  for  both 
axis.  The  area  under  each  portion  of  the  spectral  density  curve  is  the 
contribution  of  that  corresponding  frequency  to  the  variance,  and  the 
total  area  is  equal  to  the  total  variance.  Very  often  in  practice,  how¬ 
ever,  a  woiker  will  be  studying  a  variable  whose  spectrum  covers  a  range 
of  several  orders  of  magnitude  on  both  the  frequency  and  spectral  density 
axis.  In  such  a  case,  a  linear  plot  of  the  spectrum  can  obscure  much  of 
the  detail,  so  a  log-log  plot  is  often  used.  Although  the  area  under  the 
curve  is  no  longer  equal  xo  the  total  variance,  a  log-log  plot  has  the 
advantage  of  covering  wide  ranges  and  of  presenting  as  a  straight  line  any 
spectrum  which  obeys  a  power  law,  such  as  isotropic  turbulence  in  a 
stream  of  air.  A  tnird  method  of  presenting  spectra  and  a  method  which 
is  often  used,  particularly  in  meteorology,  is  to  plot  f  s(f)  against 
log  f.  Since 

Area  =  /  s(f)  df  =  /f  s(f)  d(log  f) 

the  area  under  the  curve  is  preserved,  and  one  can  still  observe  the 
relative  contribution  of  different  frequency  bands  to  the  total  variance. 
Wide  ranges  of  spectral  density  and  frequency  can  also  be  presented.  In 
Figure  A. 2. 7,  the  three  methods  of  presenting  spectra  are  illustrated  for 
a  spectrum  of  air  pressure  at  the  ground  surface  obtained  by  the  author. 


4.  Relation  of  spectral  density  function  to  autocovariance  function 
The  autocovariance  function,  R(xu),  is  defined  by 


R(t  )  = - L_ 

u 

N-l-u 


M-l-u 

2 


.I  _(H-l-u)  tp(V  '  51  ‘P(tJ  +  Tu>  '  51  A-2': 


where  the  bar  denotes  the  mean.  As  the  name  implies,  the  autocovariance 

function  expresses  a  covariance  of  p(tj)  with  itself.  The  tu  represents 

a  lag  in  time,  and  usually  it  is  an  integer  multiple  of  At,  i.e.,  tu  =  uAt 

where  u  is  an  integer.  Thus,  R(xu)  is  the  covariance  of  p(tj)  with  itself 

tu  units  of  time  later.  It  is  computed  by  multiplying  (after  subtraction 

of  the  mean)  each  observation  of  p(t.)  by  another  observation  of  p(t  ) 

J  J 

taken  a  time  x  later  (or  earlier  for  negative  xu) ,  then  summing  all  the 
lagged  products  and  dividing  by  the  number  of  products.  For  xu  ■  0, 

R(xu)  is  identical  with  the  usual  variance  function  defined  by  Equation 
A.2.16,  and  very  often  the  autocovariance  function  is  normalized  by 
dividing  it  by  the  usual  variance.  The  normalized  autocovariance  function 
is  called  the  autocorrelation  function  which  varies  only  between  ±1.  A 
plot  of  autocovariance  for  the  curve  of  Figure  A. 2. 5  is  shown  in  Figure 
A. 2. 8.  The  autocovariance  is  seen  to  be  large  and  positive  at  u  =  0 
when  the  curve  agrees  perfectly  with  itself,  small  at  u  =  2  when  the  curve 
does  not  correspond  with  itself,  large  and  negative  at  u  =  3  when  the 
positive  peaks  of  the  curve  are  opposite  the  negative  peaks ,  and  large 
and  positive  at  u  =  6  when  all  peaks  correspond.  When  N»u  and  p  =  0  (the 
data  can  be  adjusted  to  give  a  zero  mean,  if  necessary).  Equation  A. 2. 35 


becomes 


R(xu)  = 


J  =  - 


2  p(tj)  p(tj  +  xu) 

_  (»-l) 

2 


A. 2. 3 6 


ll*9 


l 


1 


\  fr 


!  I 


■  i 
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Equation  A. 2. 36  will  be  used  shortly. 

The  autocovariance  function  is  useful  because  it  illustrates  the 
lengths  of  time  for  which  a  variable  can  be  expected  to  be  correlated 
with  itself.  It  is  useful  also  because  it  happens  to  be  the  Fourier 
transform  of  the  spectral  density  function.  A  rigorous  proof  of  this 
fact,  as  given  by  Blackman  and  Tukey  (1958,  p.  72),  is  possible  only 
for  continuous  variables  extending  over  infinite  ranges.  However,  a 
demonstration  of  the  basic  idea  can  be  given  by  using  the  discrete 
equations  already  developed.  First,  consider  the  concept  of  discrete 
convolution  ,  which  is  somewhat  analogous  to  the  convolution  of  contin¬ 
uous  functions  described  by  Blackman  and  Tukey.  Let 


N-l 

o 


S(fn)  =  At 


p(t^)  exp(-i2irfntj) 


A.2.37 


N-l 

2 

and  p(tj)  =  Af  ^ 
n  =  -- 


(N-l) 


S(t  )  exp(i2irf  t . ) 
J  n  J 


A. 2. 38 


describe  a  Fourier  transform  pair.  Now  let  ( f*n )  =  S(-fn)  *  S(fn),  and 
let  the  Fourier  transform  of  ( fn )  be  denoted  B(tu). 

By  Equation  A. 2. 38 

B(tu)  =  Af  I  S^fJ  eyp(i2TTfnxu)  =  Af  T.  S(-fn)  S(fQ)  exp(i2irfnxu) 


A. 2. 39 


Then,  using  Equation  A.2.37  to  substitute  for  S(-fn) 


B(tu)  =  Af  I  [At  I  p(t.)  exp(i2irf  t.)]  S(f  )  exp(i2trf  r,)  A.2.^0 

n  i  J  njn  nu 


I 


150 


By  reversing  the  order  of  summation,  one  obtains 


B(tu)  =  At  I  [Af  E  S(fn)  exp(i2nfn(tj  +  tu))]  p(tj)  A.2.1H 


Using  Equation  A. 2. 38,  a  substitution  may  be  made  for  the  expression 
in  brackets,  so  that  Equation  A.2.41  becomes 

N-l 

2 

BCtu)  =  At  p(t  )  p(tj  +  t  )  A. 2. 42 

4  (N-l) 

J  -  -  2 


Recalling  that  T  =  (N-l)  At,  R(tu)  in  Equation  A. 2. 36  is  seen  to  be 

equal  to  the  expression  for  B(t  )  given  in  Equation  A. 2. 42  divided  by  T. 

Now  recall  that  B(tu)  was  defined  in  Equation  A. 2. 39  as  the  Fourier 

transform  of  S(-f  )  S(f  ).  Since  S(-f„)  is  the  complex  conjugate  of 
n  n  n 

.  -  rj 

S(fn),  B(tu)  is  the  Fourier  transform  of  |S(fn)r.  Using  Equation  A. 2. 21 
for  C0  =  p  =  0,  the  relation  between  the  spectral  density  function  and 

|s(fn)|2  is 

s(f  )  =  Af|s(fn)|2,  n  =  -  I2lil,-(»L3I  ,  .  .  .  ,  -1,  0,  1,  .  .  .  N-l 

2  2  ~ 

A. 2. 43 

Therefore,  using  Equation  A. 2. 43  and  the  fact  that  R(tu)  =  Af 

B( t  ) ,  one  obtains 
u 

N-l 

2 

R(t  )  =  Af  ^  s(f  )  exp(i2irf  t  )  A. 2.44 

u  Z  n  n  u 

2 

which  states  that  the  autocovariance  function  is  the  Fourier  transform  of 
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the  spectral  density  function.  Due  to  the  reciprocal  nature  of  the  Fourier 
transform  pairs,  s(fn)  must  also  be  the  Fourier  transform  of  R(tu),  i.e., 

N-l 

2 

s(fQ)  =  At  2  (N_3  \  R(tu)  exp(-i2ufnTu)  A.2.45 

u  =  -  - - — 

2 

During  the  past,  the  standard  procedure  for  computing  the  spectral 
density  function  for  a  set  of  data  points  involved  calculating  the  auto¬ 
covariance  functions  first.  Using  Equation  A. 2. 35  values  of  the  auto¬ 
covariance  function  were  computed  for  lags  covering  about  one-tenth  of  the 
total  observation  period.  The  spectral  density  function  was  then  computed 
from  Equation  A. 2. 45  for  the  number  of  lags  available.  The  development 
of  the  fast  Fourier  transform  technique  by  Cooley  and  Tukey  (1965 )  has 
made  another  approach  computationally  advantageous.  First,  the  S(fn)  are 
computed  directly  from  the  data  using  Equation  A. 2. 30  for  n  =  1,  2,  .  .  .  . 
N-l  .  Each  S(f  )  is  then  multiplied  by  twice  its  complex  conjugate  and 

-5- 

the  spectral  density  function  is  obtained  from  Equation  A. 2. 34. 

Also,  as  explained  by  Stockham  (1966),  the  autocovariance  function 
now  may  be  obtained  most  easily  not  from  lagged  products  but  by  a  second 
fast  Fourier  transformation.  After  the  spectral  density  function  is  ob¬ 
tained,  the  transformation  given  by  Equation  A. 2. 41  is  used  to  compute  the 
auto covariance  function  from  the  spectral  density  function.  However,  it 
is  not  the  usual  autocovariance  function  but  a  circular  autocovariance 
whereby  overlapping  values  at  the  end  of  the  summation  interval  for  the 
lagged  products  are  shifted  around  to  the  other  end.  This  is  why  u  goes 

from  -  (N-l)  to  (N-l)  in  Equations  A. 2. 44  and  A.2.45. 

2  2 
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Since  from  Equation  A. 2. 43,  s(fn)  can  be  seen  to  be  an  even  function 
of  f  ,  and  since  from  Equation  A. 2. 36,  R(TU)  can  be  seen  to  be  an  even 
function  of  tu,  the  summation  in  A. 2. 1+1+  and  A. 2. 1*5  need  be  carried  only 
over  half  the  ranges.  Thus, 

N-l 

2 

R(tu)  =  Af  [s ( f0 )  +  2  5  s(fn)  exP'i2irfnTu^  A.2.U6 

n=l 

and 

N-l 

2 

s(f  )  =  At  [R(t0)  +  2  5  R(t  )  exp(-i2irf  t  )]  A. 2. 1+7 

n  ^1  u  n  U 

5.  Variability  of  the  spectral  estimates 

The  experience  of  many  investigators ,  particularly  with  meteor¬ 
ological  data,  has  shown  that  the  spectral  density  estimates  calculated 
from  Equation  A. 2.3*+  may  scatter  widely  and  also  that  the  spectra  from 
two  different  portions  of  the  same  time  series  may  differ.  These  anomalies 
are  now  regarded  to  be  due  mostly  to  sampling  error,  and  methods  are 
available  to  obtain  the  underlying  "smooth"  spectrum.  The  traditional 
method  described  by  Blackman  and  Tukey  (1958)  has  been  to  calculate  the 
autocovariance  funcu  n  for  lags  extending  over  approximately  one-tenth 
of  the  total  sampling  period  and  then  obtain  the  spectrum  from  a  Fourier 
transform  of  the  autocovariance  function.  Since  the  lags  extend  only 
over  one-tenth  of  the  total  period,  the  computation  is  essentially  an 
averaging  and  smoothing  process,  and  a  smoother  spectrum  is  obtained. 
Additional  smoothing  is  obtained  by  forming  new  estimates  from  weighted 
averages,  i.e.  s;(l'n)  =  .25  s(fn_^)  +  0.50  s(fn)  +  0.25  s(fn+q)  shows 
one  set  of  weights  commonly  used.  These  weights  are  called  the  "Hanning" 


weights  after  the  man  who  first  used  them.  Ihe  method  is  also  adapted  for 
calculating  the  statistical  significance  of  the  spectral  estimates. 

Now  that  the  fast  Fourier  transform  technique  of  Cooley  and  Tukey 
(1965)  has  made  it  computationally  advantageous  to  compute  spectra 
directly  without  an  intermediate  calculation  of  an  autocovariance  function, 
another  method  of  smoothing  is  desirable.  Welch  (1967)  suggests  obtain¬ 
ing  spectral  estimates  from  several  intervals  and  averaging.  He  describes 
multiplying  the  data  by  weights  before  obtaining  the  spectral  estimates 
rather  than  multiplying  the  spectral  estimates  by  smoothing  wei^its  later. 

The  method  is  applicable  to  the  fast  Fourier  transform  techniques,  and 
like  the  autocovariance  method,  it  permits  calculation  of  the  statistical 
significance  of  the  spectral  estimates. 

If  the  data  contains  a  strong  repetitive  cycle,  such  as  the  whale 
calls  studied  by  Singleton  and  Poulter  (1967),  the  spectrum  will  contain 
a  sharp  peak,  and  another  problem  can  arise.  If  the  frequency  of  the  repeti¬ 
tive  cycle  falls  between  two  of  the  calculated  points  for  the  spectrum, 
these  two  points  and  all  the  rest  of  the  calculated  points  will  be 
affected.  As  explained  by  Bingham  et  al.  (1967),  the  effects  of  the  peak 
have  alternating  signs  and  decay  slowly  (like  — r  )  as  f  recedes  from 

I f“-nl  n 

the  peak  frequency,  f.  Before  computation  of  the  spectral  density  from 
Equation  A. 2.3^,  they  recommend  that  the  Fourier  coefficients  be  hanned 
according  to 

S'(fn)  =  0.25S(fn_1)  +  0.50S(fn)  +  0.25S(fn+1) 

when  the  time  origin  is  in  the  center  of  the  data  or  according  to 

S'(f  )  =  -0.25S(f  .  )  +  0 . 50S ( f  )  -  0.25S(f  ) 

n  n-1  n  n+1 


15* 


when  the  tiiae  origin  is  at  the  first  data  point.  The  hanning  causes  the 

effect  of  a  peak  to  decay  like  - i -  so  that  the  resulting  calculated 

Mnl3 

spectrum  will  have  a  sharper  peak  more  like  the  true  underlying  spectrum. 

An  alternate  procedure  is  discussed  by  Blackman  and  Tukey  ( 1959 ) • 

In  cases  like  that  discussed  in  the  preceding  paragraph,  where  a  strong 
repetitive  cycle  causes  a  rapid  change  of  spectral  density  with  frequency, 
they  suggest  that  the  data  be  adjusted  prior  to  the  computation  of  the 
spectrum  in  a  way  which  will  make  the  spectrum  more  flat.  Such  an 
adjustment  is  called  prewhitening  because  the  spectrum  is  made  to  resemble 
more  closely  the  flat  spectrum  of  "white"  noise  which  has  equal  spectral 
density  for  all  frequencies.  The  prewhitening  adjustment  may  be  accom¬ 
plished  by 

p'(tk)  =  p(tk)  -  a  pftj^) 

Where 


and 


p 1 ( t^ )  =  tht  adjusted  value  of  p 

a  =  a  constant  <  1.  (Blackman  and  Tukey  use  a  =  0.6 


in  one  of  their  examples.)  However,  prewhitening  by  this  means  introduces 


a  phase  change  in  the  data  which  may  be  undesirable  if  any  cross  spectral 
analysis  (soon  to  be  discussed)  is  performed. 


6.  Aliasing 

If  frequencies  higher  than  — i — were  originally  present  in  p,  they 

(2  At) 

will  unfortunately  "alias"  the  digital  data  when  it  is  sampled  at  discrete 

points,  as  illustrated  in  Figure  A. 2.9  and  explained  by  Blackman  and  Tukey 

(1958).  The  solid  line  in  the  figure  denotes  a  harmonic  whose  frequency  is 

higher  than  — i —  .  When  the  sampling  is  performed  at  the  At  intervals, 
(2At) 


gure  A. 2. 9  Illustration  of  how  digital  sampling  may  cause  a  high 
frequency  component  to  alias  and  appear  as  a  low 
frequency  component 
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tne  vari anility  in  p  due  to  this  higher  frequency  component  appears  in¬ 
stead  to  be  the  lower  frequency  harmonic  shown  by  the  dashed  line.  Hence, 
the  lower  frequency  harmonics  will  appear  larger  than  they  are  in  reality 
if  some  of  the  variability  in  p  is  due  to  harmonics  of  frequency  higier 
than  L—  .  In  the  practical  situation,  then,  one  must  make  At  small 
enough  so  that  all  important  frequencies  can  be  considered,  and  must 
filter  out  any  higher  frequency  "noise"  before  the  data  is  sampled. 
Filtering  can  be  accomplished  by  using  special  electrical  filters  ahead 
of  the  recording  equipment  or  by  using  recording  equipment  which  has  a 
long  mechanical  time  constant. 


7.  Trends 

If  the  data  show  a  general  tendency  to  increase  or  decrease  over 
the  total  sampling  period,  T,  they  are  said  to  contain  a  trend.  The 
presence  of  a  trend  means  that  a  portion  of  the  variance  is  due  to  fre¬ 
quencies  lower  than  the  lowest  di s tin guish able  frequency,  l/T.  When  the 
spectrum  is  computed,  the  spectral  density  at  the  low  frequency  end  will 
be  artificially  high  with  the  point  for  frequency  l/T  being  affected  most. 
The  usual  method  for  correction  is  to  remove  the  trend  from  the  data 
before  the  calculation  of  the  spectrum.  The  equation  of  the  trend  line 
can  be  computed  from  the  data  by  the  method  of  least  squares ;  then  the 
corresponding  trend  values  can  be  subtracted  from  each  data  point. 
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III.  CROSS-SPECTRAL  ANALYSIS 

Spectral  analysis  may  also  be  used  to  study  phase  and  amplitude 
relationships  between  the  fluctuations  of  two  entities,  p^  and  p2. 

Consider  the  Fourier  transforms,  as  defined  by  Equation  A. 2. 30,  of 
the  data  from  two  time  series.  Presume  also  that  the  data  have  zero  mean 
or  have  had  the  mean  subtracted  from  each  of  the  data  points.  Thus, 


N-l 

si  <-V  -  “  1 


j.-fcii 


p  (t  )  exp(i2mf  t  ) 
1  j  n  j 


A.  3.1 


N-l 

s2  (fn  )  =  2 

j  =  - 


(N-l) 


P2(tj)  exp(-i2TTfntj) 


A. 3.2 


_  -  (N-l)  i  n  i  N-l 

n  -  -  »»*••%  "J.*  w *  -L •  •  •  •  • 

2  2 


If  one  computes  Af  S  (-fn)  S2  (f^)  analogous  to  Equation  A. 2. 34, 
a  spectral  density  function  is  obtained,  but  it  is  quite  different  from 
the  spectral  density  function  of  Equation  A. 2. 3b.  Because  S2(fn)  is  not 
the  complex  conjugate  of  S^(-fn),  the  new  spectral  density  function  has 
both  real  and  imaginary  parts.  It  is  called  the  cross  spectral  density 
function  and  will  here  be  denoted  sc(fn). 


Jc(fn)  =  if  S1(-fn)  s2  (fn) 


A.  3.3 


That  sc(fn)  has  both  real  and  imaginary  parts  may  be  illustrated 


by  the  following.  Let 


Si(-fn)  -  +  ibn 


A. 3.4 
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S  (f  )  =  c  -  id 
2  n  n  n 


A.  3. 5 


where  a^,  Cjj,  and  b  ,  c^,  are  real  and  imaginary  parts,  respectively,  of 
S^(fn)  and  S2(fn).  Then,  using  Equations  A. 3.^  and  A. 3. 5  in  A. 3. 3, 


s-(f  )  =  f(a  c  +  b  d  )  +  i  f(b  c  -ad) 
cn  nn  nn  nn  nn 


A. 3.6 


In  general,  bncn  t  s^dg,  because  the  S1(fn)  and  Sg(fn)  are  from  two  diff¬ 
erent  time  series,  and  so  sc(fn)  is  complex. 

A  complex  number  may  also  be  written  in  the  form  of  an  amplitude 


times  a  phase  factor,  so  that 


Si(-fn)  =  A^expfi^) 
S2(fn>  -  A  exp(-i*  ) 


A. 3. 7 


A. 3.8 


where 


Aln  /an+bn 


<j>,  =  arctan  _H 


Apn  =  •'c2  +  b2 
n  n 

*n 

^2n  =  arctan  r 

These  forms  will  soon  be  useful.  Equation  A. 3.6  may  be  written 


MV  ■  C(V  -  iQ(V 


A.3.9 


The  real  part,  <2(0,  is  called  the  co-spectrum.  The  imaginary  part, 
Q(fn),  is  called  the  quadrature  spectrum.  If,  in  Equation  A. 3-6,  bncn  = 
then  »  and  the  two  times  series  are  in-phase .  Also,  then 


159 


Q(f  )  =  0,  so  s  (f  )  =  C(f  ),  and  thus,  C(f  )  is  also  called  the  in-phase 
n  c  n  n  n  - * - 

spectrum.  Recall  from  Equation  A. 2. 26  that  s(fQ)  represents  the  vaiiance 
in  p  per  frequency  increment.  Analogously,  C(f^)  represents  the  cross¬ 
covariance  between  p^  and  per  frequency  increment  for  various  fre¬ 
quencies  when  the  two  series  are  in-phase.  If,  in  Equation  A. 3.6  a&cn  = 
-bnd&,  then  ♦  =  $  +  tt/2  and  the  two  time  series  are  out-of-phase . 

Then  C(f  )  =  0,  so  sc(fn)  =  -i  Q(fn),  and  Q(fQ)  is  also  called  the  out-of¬ 
phase  spectrum.  Q(f^)  represents  the  cross-covariance  between  p^  and  p^ 
per  frequency  increment  for  various  frequencies  when  one  of  the  series 
is  shifted  exactly  1/U  period  with  respect  to  the  other. 

Hie  degree  of  phase  shift  between  the  two  time  series  for  various 
frequencies  is  measured  by  the  phase  spectrum,  defined  by 

<J>(fn)  =  arctan  A.  3.10 

C<fn) 

To  measure  the  degree  of  similarity  of  amplitude  between  the  two  time 
series  for  various  frequencies,  a  cross-amplitude  spectrum  is  defined  by 

A(fn)  =  /[C(fn)F  +  (QUn)]*  A. 3.11 

A(fn)  is  usually  normalized  and  used  in  the  form 

C(f  )2  +  Q(f  )2 

Coh  (f  )  =  - - - - —  A.  3.12 

11  sl<fn>  s2  <fn> 

Coh(fn)  is  called  the  coherence  and  measures  how  well  the  two  time  series 
are  correlated  for  various  frequencies.  As  explained  by  Panofsky  and 
Brier  (1965),  the  coherence  varies  from  0  to  1  and  is  analogous  to  the 
square  of  a  correlation  coefficient.  Probability  tables  have  been 


. 

' 


l6o 


established  for  coherence  functions  which  have  been  computed  from  trans¬ 
forms  of  cross  covariance  function,  soon  to  be  defined.  Ihe  tables  can  be 
used  to  test  the  closeness  of  the  relationship  between  the  two  series 
at  the  various  frequencies. 

However,  if  one  writes  out  Coh(fn)  explicitly  in  terms  of  raw 
Fourier  coefficients,  one  obtains 

[artVn  *  Vn>l2  ♦  E"<Vn  '  VU'2 

Coh(fn)  -  - - =  1 

(Af)  (a  -  ibj.Xajj  +  ibn)(Af)(cn  +  idn)(cR  -  i<^) 

A. 3.13 

for  every  frequency  fQ.  Although  Equation  A. 3. 13  may  seem  rather  surpris¬ 
ing  after  the  comments  in  the  preceding  paragraph,  the  explanation  is 
rather  simple.  When  computed  from  raw  coefficients.  Equation  A. 3. 12  is 
the  same  as  a  correlation  coefficient  computed  from  one  observation  pair. 
Therefore,  the  coherence  must  be  computed  from  coefficients  which  have 
been  smoothed  or  averaged  in  some  way.  The  current  practice  is  to  form  a 
weighted  average  over  several  frequencies  of  each  quantity  in  Equation 
A. 3. 12,  but,  as  discussed  by  Tick  (1967),  the  best  method  to  do  the 
smoothing  and  its  corresponding  table  of  confidence  limits  have  not  yet 
been  satisfactorily  established. 

In  Fig.  A. 3.1  are  plotted  the  coherence  and  the  phase  spectrum  for 
air  pressure  between  a  point  on  the  ground  surface  in  a  field  and  another 
point  downwind  from  the  first.  For  low  frequencies,  the  pressure  is  the 
same  at  the  two  points,  so  the  coherence  is  close  to  one  and  the  phase 
angle  between  them  is  close  to  zero.  As  frequency  increases,  the  correspond¬ 
ing  wave  lengths  of  pressure  waves  moving  across  the  field  become  smaller. 

The  pressure  at  the  upwind  point  changes  before  that  at  the  downwind  point 
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Figure  A. 3.1  Coherence  and  phase  spectrum  for  air  pressure  between 
a  point  on  the  ground  surface  in  a  field  and  another 
point  downwind  from. the  first 
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and  is  shown  in  the  figure  by  an  increasing  phase  angle  between  the  points 

with  increasing  frequency.  As  frequency  becomes  higher  (and  wave  length 

smaller),  the  two  points  do  not  always  see  the  same  wave,  as  evidenced 

by  the  decrease  in  coherence  with  frequency.  When  the  coherence  gets 

small,  the  phase  angle  changes  wildly  and  is  no  longer  meaningful. 

If  the  convolution  procedure  is  performed  on  S, (-f  )S0(f  )  exactly 

1  n  2  n 

as  it  was  on  S(-fn)S(fn)  in  Equations  A. 2. 39  through  A. 2. 42,  one  finds  that 


VTu>  ~  2  W  p2(tj  +  Tu>- 

,  _  (N-l) 


N-l 

2 


A. 3.14 


Rc  (Tu)  =  Af 


N-l 

2 

2 

n  =  ^ii 


sc(f,n)  exp(i2TifnTu)  , 


A. 3. 15 


and 


s  (f  )  =  At 
c  n 


N-l 

2 

2 

u  =  -(1,-1) 


R  ( x  )  exp(-i2irf  x  ) 
c  u  n  u 


A. 3.16 


where  R  (x..)  is  the  cross-covariance  function  and  x  is  a  lag  in  time  by 
c  u  -  u 

which  one  time  series  is  shifted  with  respect  to  the  other.  When  x^  =  0, 
Rc(xu)  is  identically  equal  to  the  usual  covariance  from  elementary- 
statistics,  i.e. 

Bc(0)  =  — ~  E  ?2  ^  ^  A. 3. IT 

J 

If  one  describes  both  p^  and  ^  Equation  A. 3. IT  by  a  Fourier  series, 
one  obtains  an  equation  analogous  to  Equation  A. 2. IT.  Proceeding  simi¬ 
larly  to  the  steps  in  Section  A. II. 2  and  Section  A.II.3*  one  finds  that 
the  contributions  of  the  various  harmonics  to  the  covariance  may  be 
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evaluated  and  plotted  to  form  e  spectrum.  The  spectrum  is  identical  to 
the  cospectrum,  and  the  area  under  the  curve  is  equal  to  the  covariance. 
The  latter  fact  can  he  used  to  evaluate  the  covariance,  and  has  been  used, 
for  instance,  by  McBean  (1968)  to  obtain  the  covariance  between  tempera¬ 
ture  and  vertical  wind  (in  this  case,  the  covariance  was  particularly 
useful  because  it  is  equal  to  the  vertical  eddy  flux  of  heat.) 

If  Equation  A. 3. 15  is  used  for  its  computation,  the  cross-covariance 
function  obtained  will  be  a  circular  cross-covariance  whereby  the  over¬ 
lapping  values  at  the  ends  of  the  time  series  due  to  the  offset  lag  iu 
are  shifted  around  to  the  front.  If  Equation  A. 3.1*1  is  used  for  it3 
computation,  the  suninring  and  averaging  are  usually  performed  only  over 
pairs  of  values  which  do  not  overlap.  One  notes  from  Equations  A. 3.15 
and  A. 3.16  that  the  cross  spectral  density  function  and  the  cross-covari¬ 
ance  function  are  Fourier  transforms  of  one  another. 


V 


1614 


IV.  TWO-DIMENSIONAL  SPECTRAL  ANALYSIS 

Some  entities  fluctuate  with  respect  to  more  variables  than  just 
time.  In  this  section,  the  concepts  introduced  for  one-dimensional 
spectral  analysis  are  expanded  to  two  dimensions.  Although  the  orientation 
taken  here  is  for  one  of  the  dimensions  to  be  time  and  for  the  other  to 
be  distance  parallel  to  the  wind,  the  concepts  apply  as  well  to  any  two 
variables.  Just  as  it  will  be  shown  how  to  extend  the  concepts  from  one 
to  two  variables,  the  concepts  can  also  be  expanded  to  include  more  than 
two  variables.  Lumley  and  Panofsky  (1961+),  for  example,  discuss  spectra 
of  entities  that  fluctuate  with  four  variables,  time  and  space  in  three 
directions. 


1.  Fourier  series  representation 

Suppose  p  represents  a  variable  such  as  pressure  in  some  turbulent 

field  moving  horizontally  over  the  ground  surface,  so  that  p  =  p(t,  x) 

where  t  it  time  and  x  is  the  horizontal  distance.  One  observes  the 

pressure  N  times  at  each  of  M  positions  in  the  turbulent  field.  At 

each  position  the  observations  are  taken  starting  at  t  =  -  |  for  a 

total  period  of  time  T.  The  observations  are  spaced  At  units  apart,  and 

the  time  at  which  each  is  taken  is  t.  =  jAt  for  j  =  ^  —  ,  .  .  .  . 

J  2 

1,  0,  1,  .  .  .  .  .  We  note  that  T  =  (N-l)At.  The  observation  posi¬ 

tions  are  spaced  along  a  line  parallel  to  the  direction  of  movement  of 


the  field.  The  total  length  covered  is  L,  and  the  first  position  is  at 


x=-—.  The  positions  are  spaced  x  units  apart,  and  an  individual 

position  is  denoted  by  x  =  £Ax  for  £  =  —  —  ,  .  .  .  .  ,-1,0,1, 

l  2 

....  “L  .  We  note  that  L  =  (M-l)Ax.  The  notation  is  illustrated  in 
Figure  A.h.l  where  the  vertical  direction  out  of  the  paper  represents  the 


magnitude  of  p. 
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Analogous  to  the  cne-dimensional  case,  p  can  be  represented  by  a 
Fourier  series  of  complex  exponentials.  Hie  two-dimensional  forms  of 
Equations  A. 2.1,  A. 2. 2,  and  A. 2. 3  are: 


p(t,  x)  =  Z 


Z  Z  C  exp[i2Tr(Si  -  2*)] 

n  m  t 

n  -  -oo  m  =  -oo  5  i  ij 


A. 4.1 


T/2  L/2 

/  /  exp[i2iT(Ei  -  S*)  ]  exp[i2v(££.  -  ^)  ]  dt  =  TL  when  q  =  -n 

-T/2  -L/2  T  L  T  L 

and  r  =  -m 


for  all  other  cases  A. 4. 2 


T/2  L/2 

c  =  —  /  /  P(t,x)  exp[i27r(^-  -  )]  dt 

n,m  TL  -T/2  -L/2  T  L 


A. 4. 3 


And,  analogously  to  the  one -dimensional  case.  Equations  A. 4.1  and 
A. 4. 3  are  discretized  to  conform  to  the  discrete  data.  Equation  A. 4.1 


becomes 

N-l 

2  2 

p(tJ5  x£)  =  2  2 

(N-l) 

n  =  - _  m  =  - 

2 


y  C  exp[i2^-^)]  A.4.4 

n  ,m  T  L 

(M-l) 


Equation  A. 4. 3  becomes 


-n,m  "  TL  -d  Z. 

.  _  (N^li  .  _  iNril 

0  2  *  2 


N-l 

2  p(t,,x£)  exp[-i2n(— -  - )]  A. 4. 5 


T  L 


As  in  the  one  dimensional  case,  the  terms  corresponding  to  n  and  -n 
represent  har-nonics  or  frequencies.  However,  now  we  have  the  additional 
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complication  of  an  added  dimension;  the  m,  -m  pairs  represent  cycles  per 
total  length  L  or  wave  numbers . 


2.  Relation  of  variance  to  spectral  density  function 
The  variance  for  p,  which  is  a  function  of  two  variables,  may  be 
written  as 


2 

s 


1 

MR-1 


N-l 

2 


2 

=  (H-l? 


N-l 

2  .  9 
2  *0  -p  J 

2 


A.  4.6 


where  the  p  is  the  mean  of  all  NM  observations. 

Note  in  Equation  A. 4.4  that  the  mean  CQ  0  assuming  MN  »  1,  sub- 
stitution  of  Equation  A. 4. 4  into  A. 4.6  yields 

s2  *  E  E  [  T.  E  C  exp[i2Tr(5li  -  lUfi  )  ]  -  C  J2  A.4.7 


j  £  n  m 


n  ,m 


o  o ' 

» 


Equation  A.4.7  is  the  numerical  integration  form  for 
T/2 


L/2 

=  —  /  /  [EEC 

TL  -T/2  -1/2  n  m 


,2  _  1 


n^expIiS^-f)] 


C0  G]  dt  dx 

J 

A. 4. 8 


Analogous  to  the  one-dimensional  case,  if  the  long  series  above  is 

written  out,  squared,  and  integrated,  all  terms  will  equal  zero  except 

those  whose  coefficients  are  C  C  or  C  C  Thus,  carry- 

-n,-m  n ,m  -n,m  n,-m 

ing  out  the  above  operations  with  the  aid  of  Equation  A. 4.2 


2 

s 


n  =  -(N-l) 
2 


M-l 

2 

2 

m  =  -(M-l) 
2 


C 

n  ,m 


,-m 


Co 


2 


o 

» 
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N-l 

2 


|C 

n  ,m 


A.  It. 9 


where  the  vertical  bars  denote  modulus.  The  squared  modulus  results 

because,  analogous  to  the  one-dimensional  case,  the  C  and  C  are 

-n  ,-m  n ,m 

complex  conjugates. 

The  term  corresponding  to  each  n,m  pair  represents  the  contribution 
of  the  n/T  th  frequency  and  m/L  th  wave  number  to  the  total  variance. 

If  we  now  define 


f  =  n/T  A. It. 10 

n 

Af  =  1/T  A. It. 11 

k  =  m/L  A.  It.  12 

m 

and  Ak  =  l/L  and  plot  A.h.13 


|Cn  m|^  against  fn  and  km>  we  obtain  a  Fourier  line  spectrum.  In  Figure 

A. It. 2  each  dot  represents  a  certain  value  of  |  C  I2  extending  in  a  line 

n  ,m 

up  out  of  the  paper. 

2 

If  t.ie  |  CQ  m|  are  distributed  uniformly  over  the  rectangle 
fn.,-_  to  fQ  and  km_-j_  to  km,  one  obtains  a  three- 
dimensional  histogram  where  the  height  of  the  rectangular  parallelpipeds 
can  be  obtained  from 


lcn,iJ2  =  (heieht)  (  Ak) 


A. It. lit 


or 


height  =  1 


C 


n,m 


2 


A. it. 15 


v 


AkAf 


M=1 

2L  * 

• 

• 

• 

• 

1 

• 

• 

• 

-<N-1) 


1 

(N-1) 

T 

2T 

•  • 

•  • 

•  • 

•  • 

'  2L»  • 

•  • 

Figure  A. 4.3  illustrates  the  3  dimensional  "periodogram. "  The 

height  of  the  rectangular  parallelpipeds  is  -  |c  I  and 

AkAf  n,m 

extends  up  out  of  the  paper. 

As  T  -*•  “  and  L  «  more  and  more  values  of  — i —  are  obtained,  the 

AfAk 

"periodogram"  becomes  a  continuous  surface.  The  smooth  continuous 
surface,  called  here  s(fn>km)  is  the  two-dimensional  spectral  density 
function.  It  is  shown  plotted  by  contour  lines  of  constant  spectral 
density  in  Figure  A. 4. 4.  Figure  A. 4. 4  is  a  two-dimensional  power  spectrum 
or  variance  spectrum  where 


s(fn;  = 


_1 _ 

AfAk 


C 


n,m 


2 


-(N-l) 
n  =  - 

2 


.  -1 ,  0 ,  +1 ,  . 


N-l 

1 


m  = 


-(M-l) 

2 


.  1,0,  +1 , 


M-l 

2 


A.4.16 


The  two-dimensional  spectral  density  function  represents  the  variance 

in  p  per  increment  of  frequency  and  wave  number.  The  total  volume  under 

the  surface  represents  the  total  variance.  The  graph  can  be  normalized,  if 

desired,  by  dividing  s(f  ,  k  )  by  the  total  variance  calculated  directly 

n  m 

from  Equation  A. 4. 6.  The  normalized  graph  will  always  have  the  volume 
under  the  graph  equal  to  one. 

The  point  where  n  =  m  -  0  is  equal  to  p^  and  represents  the  variance 

at  zero  wave  number  and  zero  frequency.  From  Equation  A. 4. 9,  however, 

note  that  it  is  not  included  in  the  summation  for  the  total  variance. 

Since  a  wave  with  a  positive  f  and  positive  k  and  a  wave  with  a 

n  m 


negative  f  and  a  negative  k  are  identical  waves  traveling  in  Lne  positive 
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Figure  A. 
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m 


i  *°U^  surfa=e  of  spectral  density  function, 
plotted  against  frequency  and  wave  number 
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x  direction,  the  wave  with  negative  f  and  km  may  be  ignored  if  the  ampli¬ 
tude  of  the  other  wave  is  multiplied  by  two.  Also,  since  a  wave  with  a 

positive  f  and  negative  k  and  a  wave  with  a  negative  f  and  positive  k 
n  m  n  m 

are  identical  waves  traveling  in  the  negative  x  direction,  the  wave  with 

the  positive  f  and  negative  k  may  be  ignored  if  the  amplitude  of  the 
n  m 

other  wave  is  multiplied  by  two. 

One  may  take  advantage  of  the  preceding  two  statements  by  defining 
wave  numbers  to  be  positive  quantities  only  and  by  regarding  frequencies 
to  be  both  positive  and  negative. 

The  decision  to  define  wave  numbers  as  only  positive  quantities 
causes  changes  in  Equations  A. 4. 9,  A.4.14,  A. 4. 15,  and  A.4.16. 

Equation  A. 4.9  becomes 


N-l 


2 


M  N-l 

2  2  2 

2  lc„,m!  -2  2 

m=0  n=0 


A.4.17 


and  Equation  A.4.16  becomes 


s 


^n’^m^ 


-(N-l) 

r  =  - 


-1,  0,  1, 


N-l 

2 


m  =  0 ,  1,  2, 


M-l 

2 


A.  4.18 


3.  Relation  of  two-dimensional  spectral  density  function  to  Fourier 
transform. 

In  manipulations  exactly  analogous  to  the  one-dimensional  case,  we 
can  form  a  discrete  two-dimensional  Fourier  transform  pair  by  substituting 


! 


i  ' 


,!! 
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the  expressions  from  C  given  by  Equation  A. 4. 4  into  Equation  A. 4. 4. 

n  ,m 

One  then  defines 


S(fn»km^  =  AtAx  Z  1  exp[-i2ir(t 

j  £ 


A.4.19 


and  then 

p(t,,x,,)  =  AfAk  Z  Z  S(f  ,k  )  exp[i2ir(t1f  -x5k_)  ]  A. 4. 20 

j  „  __  nm  j  n  *<  m 


n  m 


Equations  A.4.19  and  A. 4.20  form  the  transform  pair.  Analogous  to 
Equation  A. 2. 32  one  obtains 


S(f  ,k  )  =  ^n ,m 
n  m  AfAk 


A. 4.21 


and  also  since  C_n  _ffl  and  Cn  m  are  complex  conjugates, 

sCf  ,k  )  =  AfAk  | S(f  ,k  )|2 
n  m  1  n  m  1 


n  = 


m  = 


_  -(N-l) 
2 

-(M-l) 


,  .  .  .  .  ,  -2,  -1,  0,  1,  2,  .  .  .  . 

2 


A.  4.22 


•  •  ,  ~2 , -1,0,1,  2,  •  .  .  . 


M-l 


Thus,  the  spectral  density  surface  may  be  computed  by  a  double  Fourier 

—  (N— l)  l\r_l 

transform  of  the  data  for  n  =  — -  ,  -2,  -1,  0,  1,  2  .  .  .  .  — - 

2  2 

and  m  =  — — »  •  •  •  .  -2,  -1,  0,  1,  2  .  .  .  .  ,  and  then  multiply¬ 

ing  each  of  the  S(fn,km)  by  its  complex  conjugate  times  AfAk. 

The  computation  of  S(fn,km)  may  be  carried  out  using  the  fast  Fourier 
transform  described  by  Brigham  and  Morrow  (1967)*  To  illustrate,  let 
Equation  A. 4. 19  be  written  as : 

S(fn,km)  =  At  Z  [Ax  Z  pCt^xj^)  exp(i2TTX£km)]  exp(-i2rtjfn)  A. 4.23 

j  & 
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The  quantity  in  brackets,  denote  it  S(t.,k  )  is  simply  the  one- 

J  m 

dimensional  Fourier  transform  of  the  data.  It  is  the  same  as  Equation 
A. 2. 30  except  the  transformation  is  with  respect  to  wave  number  instead 
of  frequency.  Thus,  using  the  fast  Fourier  transform,  one  calculates 


S(tj,km)  =  Ax  I  p(t  .x^)  exp(i2nxAkm) 
£ 


A.4.24 


Then  applying  the  transform  a  second  time,  one  calculates, 


S(fn,kJ  =  At  I  S(tj,k  )  exp(-i2irkjfn) 

j 


A.  It. 25 


The  above  process  may  also  be  used  to  calculate  p(tj,x^)  from 
S(fn,km)  in  Equation  A. 4. 20. 

4.  Relation  of  two-dimensional  spectral  density  function  to  two- 
dimensional  covariance  function. 

The  two-dimensional  autocovariance  function  is  defined 


R(Tu’6v)  = 


N-l-u 

2 


(N-u)  (M-v)  -  1 


M-l-v 

2 

2 


j  =  -(N-l-u)  A  _  -(M-l-u) 
2  2 

[p(tj  +  xu,  Xy  +  6v)-p] 


fp(tj,xA)-p] 


A. 4. 26 


1 


where  the  bar  denotes  the  mean.  tu  represents  a  lag  in  time,  and  it  is 
usually  an  integer  multiple  of  At,  i.e.,  =  uAt  whe'e  u  is  an  integer. 

6y  represents  a  spacing  between  observation  stations,  and  it  usually  is 
an  integer  multiple  of  Ax,  i.e.  6y  =  vAx,  where  v  is  an  integer. 

When  N  »  u,  M  »  v,  and  p  =  0  (the  mean  can  be  subtracted  from  all 


points  if  necessary),  Equation  A. 4. 26  becomes 
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R(t  ,6  )  =  — 
U*  v'  NM 


N-l 

2 


2 

=  -(N-ii 


M-l 

2 

2  P(tj,  x4)  p(tj  +  tu,  x£  +  6y) 

v  =  -(M-i) 

2 


A. 4. 27 


A  convolution  procedure  analogous  to  the  one-dimensional  case  reveals 
the  two-dimensional  autocovariance  function  to  he  the  double  Fourier 
transform  of  the  two-dimensional  spectral  density  function  so  that 


R(iu,6v)  =  AfAk 

n 


N-l 

2 


M-l 

2 

2  s(fn>km*  exp[i2ir(fnTu  -  km6y)]  A.4.28 

-(M-l) 


The  spectral  density  function  must  also  be  the  double  Fourier  trans¬ 
form  of  the  autocovariance  function  or 


3(fn*km)  =  AtAx 


u  = 


-(N-l) 

2 


M-l 

2 

2  r^tu»6v)  exp[-i2Tr(fnxu  -  y>v)]  A. It. 29 

.  -(M-l) 

2 


The  autocovariance  function  in  Equation  A.  It. 28  is  not  the  usual  auto¬ 
covariance  function  but  a  circular  autocovariance  function  where  the  over¬ 
lapping  values  at  the  ends  of  the  summation  intervals  in  the  usual  auto¬ 
covariance  function,  Equation  A. It. 28, have  been  shifted  to  the  fronts. 

This  is  why  u  goes  from  ~--R~ 3:)-  to  and  v  from  — (%ii  to  HlzL  in 

2  2  2  2 

Equations  A. It. 28  and  A.  4.29  even  though  the  assumption  was  stated  that 
N  »  u  and  M  »  v  in  going  from  Equation  A. 4.26  to  Equation  A. 4. 27. 

It  can  be  deduced  from  Equation  A. 4.27  that  R(xu,5u)  =  r(-tu,-6v)  and 
that  R(-tu,6v)  =  R(tu,-6v).  Thus ,  Equation  A. 4.29  may  be  modified  to 
include  only  positive  spacings  so  that 
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s(fn»km)  =  2AtAx 


N-l 

2 


2 


=  -(N-l) 


2 


M-l 

2 

2  R(tu*6v)  exp[-i2ir(fnTu  -  k^)] 

V  =  0 


A.U.30 


Prom  Equation  A.lt.18  it  can  also  be  seen  that  Equation  A.l{.28  may 
be  modified  to  include  only  positive  wave  numbers  so  that 


R(tu,6v)  =  2AfAk 


N-l 

2 


n 


2 

_  -(N-l) 


2 


M-l 

2  s<fn'km>  exp[i2T.(rnTu  -  Vt)1 
m  =  0 
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V.  USE  OF  CROSS  SPECTRAL  DENSITY  FUNCTIONS  IN  TWO-DIMENSIONAL  ANALYSIS 
To  calculate  the  two-dimensional  spectral  density  function  from 
Equations  A.U.19  and  A.U.22,  N  observations  are  required  at  each  of  M 
observation  stations.  Since  N  and  M  are  both  usually  very  large  numbers, 
it  is  desirable  to  have  methods  which  require  less  acquisition  and  process¬ 
ing  of  data.  Under  certain  assumptions,  discussed  in  this  section,  it  is 
possible  to  calculate  the  two-dimensional  spectral  density  function  from 
data  taken  at  drastically  fewer  than  M  stations.  The  orientation  is  again 
toward  an  entity  which  fluctuates  with  time  and  one  direction  in  space, 
but  the  methods  can  be  applied  wherever  the  assumptions  are  valid. 


1.  The  assumptions  of  stationarity  and  spacial  homogeneity 
Suppose  p  is  observed  with  time  at  several  points  separated  by  a 
spacing  =  vAx  units  of  distance  from  some  space  origin.  The  sampling 
scheme  is  illustrated  in  Figure  A.5«l-  One  subtracts  the  mean  from  all 
observations  so  that  only  the  fluctuations  in  p  are  considered. 

When  p  is  statistically  stationary,  the  mean,  the  variance,  and  other 
statistical  parameters  of  p  do  not  change  with  time  even  though  p  itself 
is  fluctuating  constantly.  If  p  is  spacially  homogeneous,  the  mean,  the 
variance,  and  other  statistical  parameters  do  not  change  with  distance. 
When  one  can  assume  that  p  is  statistically  stationary  and  spacially 
homogenous ,  the  two-dimensional  autocovariance  function 


N-l  M-l 

2  2 

R(tu>6v>  =  ^  2  2  p(tj,x£)p(tj  +xu,  x£  +  6U) 

1  -  -(N-l)  .  _  -OM  ) 


A. 5.1 


is  a  function  only  of  the  time  and  space  increments  x  and  6  .  It  does 

u  v 

not  change  with  a  new  origin  of  time  or  a  new  origin  of  space. 


Figure  A. 5.1  Sampling  scheme  which  can  he  used  where  p  is  statistically 
stationary  and  spacially  homogeneous 
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Equation  A. 5.1  can  be  written 

N-l 

B(V4t>  '  M  !N  2  pCtj*X-(N-l)  1  p(tj  *  7«>’  x-(N-l)  *  6v> 

,  -(M-D  2  2 

J  -  2 

+ . 

+ . 

+ . 

N-l 

,  2 

+f  2  P(tj,x0)p(tj  +  Tu,  X0  +  5V) 

-(N-l) 

J  2 

+ . 

+ . 

+ . 

+ . 

+ . 

N-l 

+  !  2  p(vxk_1wvtu-\-i  +  6v))  a-5-2 

.  .  -(N-D..  —  ~T 

J  2 


and  one  notes  that  under  the  above  assumption,  all  of  the  sums  on  j  are 
equal.  Since  there  are  M  of  the  sums, 

N-l 

R(VSv>  =  N  2  Pttj *  v  x0  ♦  «v)  A.5.3 

_  -(N-l) 

J  "  2 

Recalling  Equation  A. 3.1*+  one  sees  that  Rtx^jS^)  is  equal  to  the  cross 
covariance,  RC(TU»<$V)>  between  the  pressure  at  point  xQ  snd  point 
xQ  +  6  .  Thus , 

°  y 

R(Vsv!  -  R«.(tu,5v) 


A.5.*+ 


which  is  an  important  simplification  and  can  drastically  reduce  the 
required  number  of  measurements  where  the  assumptions  are  valid. 


2.  Relationship  between  cross  spectral  density  function  and  2-D 
spectral  density  function. 

The  Fourier  transforms  of  the  data  gathered  by  the  scheme  in  Figure 
A. 5.1  can  be  computed  for  each  of  the  sampling  points  according  to  Equation 
A. 2. 30  which  is  rewritten  here. 


M-l 

2 

S(-f  ,x0)  =  At  ^  p(t  ,x0)  exp(i2irfnt.)  A. 5. 5 

-(M-l) 

J  "  2 

N-l 

2 

S(fn>  xQ  +  6y)  =  At  ^  p(t  ,  xD  +  5y)  exp(-i2irfntj)  A. 5. 6 

_  -(N-l) 

J  -  2 


The  cross  spectral  density  function  between  the  two  points,  recalling 
Equation  A. 3.. 3,  is 


sc(fn,6v)  =  Af  S(-fn,x0)s(f  ,  xo  +  6  )  A. 5. 7 

By  Equation  A. 3.16,  it  is  also  equal  to  the  Fourier  transform  of 

the  cross-covariance  function 

N-l 

~2 

sc(fn>6v)  =  At  2  Rc(tu,6v)  exP(i2lTfnTu)  A*5’8 

-(N-l) 

u  =  2 

Now,  suppose  one  computes  a  two-dimensional  spectrum,  called  here  W(fn,km), 

by  a  Fourier  transform  of  sc(fn,Sv)  according  to: 

M-l 
2 

W(fnskm^  =  Ax  2  sC(fn»6v)  exp(i27rkmov) 

-(M-l) 


A. 5. 9 


182 


If  one  substitutes  for  sc^n»^y)  fro®  Equation  A. 5. 8  into  A.5.9j  one  obtains 


M-l  N-l 

2  2 

^n>K'  =  I  x  At  2 

-(N-l) 


R  (iu>6v)exp(-i2nfnTu)  exp(i2irkm6v) 


v  =  2 


u  = 


A. 5. 10 


which  can  be  rewritten  to  yield 


N-l  M-l 

2  2 

=  AtAx  2  2  Rc(tu>6v)  exp[-i2.(fnTli  -yy)] 

-(N-l)  -(M-l) 

u  =  2  v  =  - 2 

A. 5. 11 

Recalling  Equations  A. 5-^  and  A.1+.29  one  can  see  that  W(fn,km)  equals 
s(fn,km).  Thus,  the  two-dimensional  spectral  density  function  may  be 
computed  by  a  double  Fourier  transform  of  the  two-dimensional  covariance 
function,  i.e., 

N-l  M-l 

2  2 

s(f  ,km)  =  AtAx  ^  2  R(Tu,6v)exp[-i2Ti(fnTu  -km6v)]  A.5.12 

All  -(M-l) 

u  =  2  ▼  •  - 1 — 


or  it  also  may  be  computed  by  a  single  transform  from  the  cross  spectral 
density  function,  i.e.. 


M-l 

2 

s(f  ,k  )  =  Ax  2 

n  m  *- 


sc(fn,6v)exp(i2^km6v) 


A. 5. 13 


-(M-l) 
v  -  — - - - 


2 
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3.  Relationship  between  cross  spectral  density  functions  and  two- 

dimensional  auto-covariance  function. 

Suppose  one  computes  a  two-dimensional  function,  called  here  Y(f  ,6v), 

by  a  Fourier  transform  of  s  (f  ,iS  )  according  to: 

c  n  v 


Y(tu,6v)  =  Af  $  sc(fn,6v)exp(i27rfnTu) 
-(N-l) 


A.5.1J+ 


The  inverse  of  Equation  A. 5. 13, 


sc(fn’6v)  =  ^  2  s(f  .kje^^ik^) 

-(M-l) 

m  = 


may  be  used  to  substitute  in  Equation  A.5.14  to  obtain: 


A.5.15 


Y(V6v}  =  Af  ,  ,  M  ,  sCfn,km)exp(i2^km6u)  exp(i27rfnTu) 

n  =  -(N-l)  m  =  ~  VW--L  i 


n  =  m  =  ZSjtl 

2  2 


A. 5-16 


Upon  rearrangement,  the  right-hand  side  is  seen  to  be  identical  to  the 
right-hand  side  of  Equation  A. 4.28.  Therefore,  Y(xu,6v)  =  R(tu,6v)  or, 
again  writing  A. 4.28, 


N-l  M-l 

2  2 

B<vV  ■ 4f4k  2  2 

n  *  -tgr.il  „  .  -(M-l) 


s(f  ,k  )exp[+i2ir(f  t  -  k  6  )]  A. 5. 17 
n  m  n  u  m  v 


Thus,  the  two-dimensional  auto- covariance  may  be  computed  from  a  double 
Fourier  transform  of  the  two-dimensional  spectral  density  function,  or  it 
may  be  computed  by  a  single  Fourier  transform  of  the  cross-spectral  density 
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function  according  to  Equation  A.5.11*,  rewritten  as: 

N-l 

2 

R(VV  =  Af  5  Sc(fn,6v)exp(i2lTfnTu)  A.5.18 

=  zlJbll 

C. 


n 
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VI.  SUMMARY 

Hie  underlying  principles  of  spectrum  analysis  are  discussed. 
Beginning  with  the  definition  of  a  Fourier  series,  it  is  shown  that 
a  spectrum  is  a  plot  of  the  squared  amplitudes  versus  frequency  of  the 
harmonics  which  contribute  to  the  variance  a  time  series.  It  is  also 
described  how  the  squared  Fourier  amplitudes  are  proportional  to  the 
spectral  density  function  computed  by  a  Fourier  transform  of  the  data. 

Hie  spectral  density  function  is  also  shown  to  be  the  Fourier  transform 
of  the  more  familiar  auto-covariance  function. 

Hie  discussion  next  considers  two  simultaneous  mime  series  and  shows 
how  the  closeness  of  the  relationship  of  and  the  phase  angle  between 
the  two  series  may  be  assessed  by  computation  of  the  cross  spectral 
density  function  from  Fourier  transforms  of  the  data  from  both  series. 

Hie  cross  spectral  density  function  is  found  to  be  the  Fourier  transform 
of  the  <_  ^variance  function  for  the  two  series. 

Two-dimensional,  spa^e-time  spectral  analysis  is  considered.  Analo¬ 
gous  to  the  one -dimensional  case,  it  is  shown  that  a  two-dimensional 
spectrum  is  a  contour  surface  of  squared  Fourier  amplitudes  versus 
frequency  and  wave  number  of  the  harmonics  contributing  to  the  variance 
of  some  entity  in  time  and  space.  The  squared  Fourier  amplitudes  are 
shown  to  be  proportional  to  a  two-dimensional  spectral  density  function 
which  may  be  computed  by  a  double  Fourier  transform  of  the  data.  Hie  two- 
dimensional  spectral  density  function  is  shown  to  be  the  Fourier  trans¬ 


form  of  a  two-dimensional  auto-covariance  function. 
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Lastly,  it  is  shown  that,  if  samples  of  time  series  are  available 
for  various  points  across  a  field,  both  the  two-dimensional  spectral 
density  function  and  the  two-dimensional  auto-covariance  function  can  be 
computed  from  the  cross  spectral  density  functional  relationships  among 
the  points. 


187 


VII.  LITERATURE  CITED 

AJlen,  L.  H. ,  Jr.  1968.  Turbulence  and  wind  speed  spectra  within  a 
Japanese  larch  plantation.  J.  of  Appl.  Met.  7:73-78. 


Bingham,  C.,  Godfrey,  M.  D. ,  and  Tukey,  J.  W.  1967.  Modem  Techniques 
of  Power  Spectrum  Estimation.  IEEE  Transactions  on  Audio  and 
Electroacoustics  AU-15(2):56-66. 


Blackman,  R.  B.  and  Tukey,  J.  W.  1958.  The  Measurement  of  Power  Spectra. 
Dover  Publications,  New  York,  190  p. 


Bri$iam,  E.  0.  and  Morrow,  R.  E.  1967.  The  Fast  Fourier  Transform.  IEEE 
Spectrum  l+(3):63-70. 

Cooley,  J.  W.  and  Tukey,  J.  W.  1965-  An  Algorithm  for  the  Machine 

Calculation  of  Complex  Fourier  Series.  Math.  Comput.  19:297-301. 

Desjardins,  R.  L.  1967.  Time  Series  Analysis  in  Agrometeorological 
Problems  with  Emphasis  on  Spectral  Analysis.  Canadian  J.  Plant 
Sci.  1*7;  1+77-1+91. 


Jenkens,  G.  M.  and  Watts,  D.  G.  1968.  Spectral  Analysis  and  Its 
Applications.  Holden-Day,  San  Francisco,  525  pp. 


Lumley,  J.  L.  and  Panofsky,  H.  A.  1961*.  The  Structure  of  Atmospheric 
Turbulence.  John  Wiley  and  Sons,  Inc.,  New  York,  239  pp. 


McBean,  G.  A.  1968.  An  investigation  of  turbulence  within  the  forest. 
J.  of  Appl.  Met.  7:1+10-1*16. 


188 


Panofsky,  H.  A.  and  Brier,  G.  W.  1958.  Some  Appli cation  of  Statistics 

to  Meteorology.  The  Pennsylvania  State  Univ. ,  Univ.  Park,  Penn,  223  p. 

Rodriguez-Iturbe,  S.  and  Nordin,  C.  F.  1969.  Some  applications  of 

cross-spectral  analysis  in  hydrology:  rainfall  and  runoff.  Water 
Resources  Res.  5:608-621. 

Singleton,  R.  C.  and  Poulter,  T.  C.  1967.  Spectral  analysis  of  the 
call  of  the  male  killer  whale.  IEEE  Transactions  on  Audio  and 
Electroacoustics  AU-15 (2 ): 10^-113. 

Stockham,  T.  G.  1966.  High-speed  convolution  and  correlation.  Spring 
Joint  Computer  Conference.  American  Federation  of  Information 
Processing  Societies  Proc.  28:229-233.  Spartan  Books,  Washington, 

D.  C. 

Tick,  L.  J.  1967-  Estimation  of  coherency.  In  Harris,  B.  ed. ,  Spectrum 
Analysis  of  Time  Series.  John  Wiley  and  Sons,  Inc.,  New  York,  319  pp. 

Welch,  P.  D.  1967.  The  use  of  the  fast  Fourier  transform  for  the 

estimation  of  power  spectra:  A  method  based  on  time  averaging  over 
short,  modified  periodograms .  IEEE  Transactions  on  Audio  and 
Electroacoustics  AU-15(2) : 70-73. 


UflCLASfiTFTF.n  - 

Security  CU«»ific*tion 


DOCUMENT  CONTROL  DATA  •  R  &  D 


I.  ORIGINATING  ACTIVITY  (Corporal*  author) 

Microclimate  Investigations,  SWC-ARS-USDA 

Bradfield  Hall;  Cornell  University 

Ithaca.  New  York _ lliflgn _ _ _ 

24.  REPORT  SECURITY  C  LA  SSI  PIC  ATION 

mrxAggjFXffl? 

26.  GROUP 

3  REPORT  TITLE 

EFFECTS  OF  AIR  TURBULENCE  UPCN  GAS  EXCHANGE  FROM  SOIL 

with  appendix 

RARTf!  rrmPV.PTR  f!P  RPPP*PPAT  avatvotc!  ov  rtATiRi r  imum 

4,  desc  ripti  ve  NOTES  (Typo  of  rapori  and  Inelualwa  data*) 

Interim  Reoort 

t.  author:* l  (Urol  noma,  mlddto  Initial,  tool  noma) 

B.  A.  Kimball  and  E.  R.  Lemon 

«.  REPORT  OATE 

November  1969 

7 A.  TOTAL  NC.  OP  PAOES  7b.  NO.  OP  REPS 

»«  CONTRACT  OR  GRANT  NO. 

Cross-Service  Order  2-63 

b.  PROJECT  NO. 

«.  DA  Task  1T0-61102-B53A-17 

d. 

M.  ORIGINATOR'*  REPORT  NUMBERIS)  ' 

USDA-ARS-SWC  No.  405  and  406 

I. U. Research  Report  No.  879  and  870A 

Ma  report) 

EC0M  2-68I-U  and  ECOM  2-681-5 

10.  OiSTRIBUTION  STATEMENT  £ 

Distribution  of  this  document  is  unlimited 

11.  SUPPLEMENTARY  NOTES 

12.  SPONSORINB  MILITARY  ACTIVITY 

U.  S.  Amy  Electronics  Command 

Atmospheric  Sciences  Laboratory 

Fort  Huaehiien  -  A-rt  7T>nn _ _ 

13.  ABSTRACT 
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soil  or  ether  porous  media  wac  constructed.  Subsequent  measured  rates  of  znr,  move¬ 
ment  from  beneath  surface  coverings  of  coarse  and  fine  gravel,  very  coarse  and 
medium  sand,  Chenango  silt  loam  and  s t raw  were  correlated  against  concurrent 
measurements  of  wind  and  air  pressure  fluctuations.  At  a  depth  of  2  cm  gas  movement 
in  silt  loam  va3  not  significantly  affected  by  air  turbulence  but  in  coarse  gravel 
was  increased  10-fold  by  wind  velocities  of  U00  cm/sec.  Intermediate  effects  were 
observed  for  media  of  intermediated  none  sizes.  1  /\' 

Mathematical  equations  for  calculation  of  soil  air  mass  flow  are  derived,  They 
are  based  upon  the  spectrum  of  air  pressure  fluctuations  at  the  ground  surface.  Air 
pressure  in  the  field  wan  measured  during  both  day  and  night  with  and  without  a 
corn  crop.  The  spectra  of  air  pressure  were  calculated  and  all  were  roughly  straight 
lines  with  a  slope  of  about  -6/3  on  a  log-log  plot.  Spectral  density  decreased  from 
10°  to  10"-Lubar2 /cycle /sec  over  a  frequency  range  from  IQ*”14  to  10s  cyclc/sec.  Soil 
air  mass  flow  was  calculated  from  the  spectra  for  several  soils.  Using  previous 
work,  an  attempt  was  made  to  evaluate  the  increase  of  soil  gas  movement  beyond 
diffusion  from  the  soil  air  mass  flow.  *1116  predicted  increases  in  soil  gas  movement 
were  lower  than  the  observed  increases;  several  reasons  fer  the  discrepancy  are 
discussed.  / ' 

Spectral  density,  Fourier  transformation,  auto-covariance,  cross  spectral 
density,  cross-covariance,  and  other  concepts  of  spectral  analysis  are  discussed 
at  an  elementary  level  in  the  appendix. 
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